/**
 * Autor:   Lukas Brzobohaty, xbrzob06@stud.fit.vutbr.cz
 *          David Marecek, xmarec10@stud.fit.vutbr.cz
 * Soubor:  main.cc
 * Projekt: Implementace modelu zalozeneho na metode Monte-Carlo
 * Predmet: IMS
 * Popis:   
 * Poznamka:
 **/

#include "simlib.h"
#include "generator.h"
#include <iostream>
#include <vector>
#include <math.h>
#include <time.h>
#include <ctype.h>
#include <algorithm>
#include <fstream>
#include <regex.h>
#include <limits.h>

using namespace std;

/**
* Funkce pro vypocet prumeru a rozptylu
*
*
**/

void vypocet()
{
	srand(time(NULL));
	
	// inicializace pocatecnich hodnot
	double H_0 = 0.84;
	double Z_0 = 0.81;
	double M = 1.5;
	double M_t = 1.35;
	double A_s = 66.672;
	double Q_in = 0.0232; 
	double Q_f = 0;
	double delta_t = 0.1;
	double t = 0;
	double H_d = 0.84;
	double B = H_d;
	double H_k = H_d;
	
	for(int i = 0; i < POCET_KONSTANT; i++)
	{
		seeds[i] = rand();
	}
	
	for(int i = 0; i < POCET_GEN_VZORKU; i++)
	{
		double b_0 = Generator(0)*(0.1-0.03)+0.03;
		double s_0 = Generator(1);
		double alfa_1 = Generator(2)*(0.01-0.008)+0.008;
		double alfa_2 = Generator(3)*(0.01-0.008)+0.008;
		double alfa_3 = Generator(4)*(0.01-0.008)+0.008;
		double beta_1 = Generator(5);
		double beta_2 = Generator(6);
		double beta_3 = Generator(7);
		
		double Q_bmax = M * (b_0 * pow(H_k, 3/2) + 0.4 * ((B - b_0)/H_d) * pow(H_k, 5/2));	// kulminacni prutok vody natrzi
	
		double Q_b = M * b_0 * pow(H_0 - Z_0,3/2) + M_t * s_0 * pow(H_0 - Z_0, 5/2);	// prutok vody natrzi
		H_0 = (delta_t / (A_s * H_0)) * (Q_in - Q_b - Q_f) + H_0;						// uroven hladiny v nadrzi

		double A_b = b_0 * (H_0 - Z_0) + s_0 * pow(H_0 - Z_0, 2);						// prutocny prurez natrzi
		double v_0 = Q_b / A_b;															// prurezova rychlost proudeni vody natrzi

		Z_0 = -delta_t * alfa_1 * pow(v_0, beta_1) + Z_0;								// nejvyssi kota dna natrze hraze

		// vypis vysledku
		Print("%.6f\t%.6f\t%.6f\t%.6f\t%.6f\t%.6f\t%.6f\t%.6f\n", t, Q_bmax, Q_b, H_0, b_0, alfa_1, alfa_2, alfa_3);
	
		if(Z_0 > 0)
		{
			b_0 = delta_t * alfa_2 * pow(v_0, beta_2) + b_0;							// sirka natrze
		}
		else if(Z_0 == 0)
		{
			b_0 = delta_t * alfa_3 * pow(v_0, beta_3) + b_0;							// sirka natrze
		}
		else 
		{
			Z_0 = 0;
		}
		
		t += delta_t;
	}

}

/**
* Hlavni funkce
*
*
**/
int main()
{
	SetOutput("ims.dat");
  	Print("# Time\t\tQbmax\t\tQb\t\t\tH0\t\t\tb0\t\t\talfa1\t\talfa2\t\talfa3\n");
  	vypocet();
  	fstream data;
	
	double Qbmax_tabmax = (double)LONG_MIN, Qbmax_distr_tabmax = (double)LONG_MIN, b_0_tabmax = (double)LONG_MIN, alfa_1_tabmax = (double)LONG_MIN, alfa_2_tabmax = (double)LONG_MIN, alfa_3_tabmax = (double)LONG_MIN;
	double Qbmax_tabmin = (double)LONG_MAX, Qbmax_distr_tabmin = (double)LONG_MAX, b_0_tabmin = (double)LONG_MAX, alfa_1_tabmin = (double)LONG_MAX, alfa_2_tabmin = (double)LONG_MAX, alfa_3_tabmin = (double)LONG_MAX;
	double Qbmax_soucet = 0, Qbmax_prum, Qbmax_rozpt = 0, b_0_soucet = 0, b_0_prum, b_0_rozpt = 0, alfa_1_soucet = 0, alfa_1_prum, alfa_1_rozpt = 0, alfa_2_soucet = 0, alfa_2_prum, alfa_2_rozpt = 0, alfa_3_soucet = 0, alfa_3_prum, alfa_3_rozpt = 0;
	
	string line;	// string pro nacteni jednoho radku ze souboru
	data.open("ims.dat", ios::in | ios::out);
	
	if(data)
	{
		while(getline(data, line))	// cteni souboru
		{
			if(line.find("#") != 0)
			{
				regmatch_t pmatch[6];
			  	regex_t re;                                        // inicializace regexp
			  	int status;
			  	
			  	char pattern[]="^.*\\s(.*)\\s.*\\s.*\\s(.*)\\s(.*)\\s(.*)\\s(.*)$"; // maska nactena data
			  
			  	if(regcomp( &re, pattern, REG_EXTENDED)!= 0)	// kompilace regexpu
			  	{
					cerr << "Error while compiling regexp !" << endl;
					return EXIT_FAILURE;
			  	}
			  
			  	status=regexec( &re, line.data(), 6, pmatch, 0); // kontrola zadaneho parametru
			  	if(status == 0)
			  	{
			  		// vypocet max. hodnot
					Qbmax_tabmax = max(Qbmax_tabmax, atof((line.substr((pmatch[1].rm_so), (pmatch[1].rm_eo - pmatch[1].rm_so))).c_str()));
					b_0_tabmax = max(b_0_tabmax, atof((line.substr((pmatch[2].rm_so), (pmatch[2].rm_eo - pmatch[2].rm_so))).c_str()));
					alfa_1_tabmax = max(alfa_1_tabmax, atof((line.substr((pmatch[3].rm_so), (pmatch[3].rm_eo - pmatch[3].rm_so))).c_str()));
					alfa_2_tabmax = max(alfa_2_tabmax, atof((line.substr((pmatch[4].rm_so), (pmatch[4].rm_eo - pmatch[4].rm_so))).c_str()));
					alfa_3_tabmax = max(alfa_3_tabmax, atof((line.substr((pmatch[5].rm_so), (pmatch[5].rm_eo - pmatch[5].rm_so))).c_str()));

					// vypocet min. hodnot
					Qbmax_tabmin = min(Qbmax_tabmin, atof((line.substr((pmatch[1].rm_so), (pmatch[1].rm_eo - pmatch[1].rm_so))).c_str()));
					b_0_tabmin = min(b_0_tabmin, atof((line.substr((pmatch[2].rm_so), (pmatch[2].rm_eo - pmatch[2].rm_so))).c_str()));
					alfa_1_tabmin = min(alfa_1_tabmin, atof((line.substr((pmatch[3].rm_so), (pmatch[3].rm_eo - pmatch[3].rm_so))).c_str()));
					alfa_2_tabmin = min(alfa_2_tabmin, atof((line.substr((pmatch[4].rm_so), (pmatch[4].rm_eo - pmatch[4].rm_so))).c_str()));
					alfa_3_tabmin = min(alfa_3_tabmin, atof((line.substr((pmatch[5].rm_so), (pmatch[5].rm_eo - pmatch[5].rm_so))).c_str()));
					
					// soucet hodnot
					Qbmax_soucet += atof((line.substr((pmatch[1].rm_so), (pmatch[1].rm_eo - pmatch[1].rm_so))).c_str());
					b_0_soucet += atof((line.substr((pmatch[2].rm_so), (pmatch[2].rm_eo - pmatch[2].rm_so))).c_str());
					alfa_1_soucet += atof((line.substr((pmatch[3].rm_so), (pmatch[3].rm_eo - pmatch[3].rm_so))).c_str());
					alfa_2_soucet += atof((line.substr((pmatch[4].rm_so), (pmatch[4].rm_eo - pmatch[4].rm_so))).c_str());
					alfa_3_soucet += atof((line.substr((pmatch[5].rm_so), (pmatch[5].rm_eo - pmatch[5].rm_so))).c_str());
			  	}
			  	regfree( &re);
			}	
		}
	}
	data.close();
		
	// vypocet vyberoveho prumeru
	Qbmax_prum = Qbmax_soucet/POCET_GEN_VZORKU;
	b_0_prum = b_0_soucet/POCET_GEN_VZORKU;
	alfa_1_prum = alfa_1_soucet/POCET_GEN_VZORKU;
	alfa_2_prum = alfa_2_soucet/POCET_GEN_VZORKU;
	alfa_3_prum = alfa_3_soucet/POCET_GEN_VZORKU;

	data.open ("ims.dat", ios::in | ios::out);

	if(data)
	{
		while(getline(data, line))	// cteni souboru
		{
			if(line.find("#") != 0)
			{
				regmatch_t pmatch[6];
			  	regex_t re;                                        // inicializace regexp
			  	int status;
			  	
			  	char pattern[]="^.*\\s(.*)\\s.*\\s.*\\s(.*)\\s(.*)\\s(.*)\\s(.*)$"; // maska nactena data
			  
			  	if(regcomp( &re, pattern, REG_EXTENDED)!= 0)	// kompilace regexpu
			  	{
					cerr << "Error while compiling regexp !" << endl;
					return EXIT_FAILURE;
			  	}
			  
			  	status=regexec( &re, line.data(), 6, pmatch, 0); // kontrola zadaneho parametru
			  	if(status == 0)
			  	{
					// rozptyl
					Qbmax_rozpt += pow(atof((line.substr((pmatch[1].rm_so), (pmatch[1].rm_eo - pmatch[1].rm_so))).c_str()) - Qbmax_prum, 2);
					b_0_rozpt += pow(atof((line.substr((pmatch[2].rm_so), (pmatch[2].rm_eo - pmatch[2].rm_so))).c_str()) - b_0_prum, 2);
					alfa_1_rozpt += pow(atof((line.substr((pmatch[3].rm_so), (pmatch[3].rm_eo - pmatch[3].rm_so))).c_str()) - alfa_1_prum, 2);
					alfa_2_rozpt += pow(atof((line.substr((pmatch[4].rm_so), (pmatch[4].rm_eo - pmatch[4].rm_so))).c_str()) - alfa_2_prum, 2);
					alfa_3_rozpt += pow(atof((line.substr((pmatch[5].rm_so), (pmatch[5].rm_eo - pmatch[5].rm_so))).c_str()) - alfa_3_prum, 2);
			  	}
			  	regfree( &re);
			}	
		}
		
		// vypocet vyberoveho rozptylu
		Qbmax_rozpt = Qbmax_rozpt/POCET_GEN_VZORKU;
		b_0_rozpt = b_0_rozpt/POCET_GEN_VZORKU;
		alfa_1_rozpt = alfa_1_rozpt/POCET_GEN_VZORKU;
		alfa_2_rozpt = alfa_2_rozpt/POCET_GEN_VZORKU;
		alfa_3_rozpt = alfa_3_rozpt/POCET_GEN_VZORKU;
		
	}
	data.close();

	SetOutput("ims_tab.dat");
	
	Print("#\t\t\tB0\t\talfa1\talfa2\talfa3\tQbmax\n");
	Print("max\t\t\t%.5f\t%.5f\t%.5f\t%.5f\t%.5f\n", b_0_tabmax, alfa_1_tabmax, alfa_2_tabmax, alfa_3_tabmax, Qbmax_tabmax);
	Print("min\t\t\t%.5f\t%.5f\t%.5f\t%.5f\t%.5f\n", b_0_tabmin, alfa_1_tabmin, alfa_2_tabmin, alfa_3_tabmin, Qbmax_tabmin);
	Print("vyb.prum\t%.5f\t%.5f\t%.5f\t%.5f\t%.5f\n", b_0_prum, alfa_1_prum, alfa_2_prum, alfa_3_prum, Qbmax_prum);
	Print("vyb.so\t\t%.5f\t%.5f\t%.5f\t%.5f\t%.5f\n", sqrt(b_0_rozpt), sqrt(alfa_1_rozpt), sqrt(alfa_2_rozpt), sqrt(alfa_3_rozpt), sqrt(Qbmax_rozpt));
}
